﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace System.Windows.Media.Media3D10
{
    /// <summary>
    /// Classe permettant de réaliser des transformée de Fourier rapide.
    /// </summary>
    public static class FastFourierTransform
    {
        #region Public Methods
        /// <summary>
        /// Réalise une transformation de Fourier 2D
        /// </summary>
        /// <param name="input">Tableau de <see cref="Complex"/> à transformer</param>
        /// <returns>Tableau de <see cref="Complex"/> transformé</returns>
        public static Complex[,] ComputeForward(Complex[,] input)
        {
            return Compute(input, false);
        }

        /// <summary>
        /// Réalise une transformation inverse de Fourier 2D
        /// </summary>
        /// <param name="input">Tableau de <see cref="Complex"/> dans le domaine de Fourier</param>
        /// <returns>Tableau de <see cref="Complex"/> transformé</returns>
        public static Complex[,] ComputeInverse(Complex[,] input)
        {
            return Compute(input, true);
        }

        /// <summary>
        /// Obtient le module d'une transformée de fourier 
        /// </summary>
        /// <param name="input">Table de <see cref="Complex"/> représentant la transformée de Fourier</param>
        /// <returns>Tableau de <see cref="Complex"/> représentant le module de la transformée</returns>
        public static double[,] GetMagnitudePlot(Complex[,] input)
        {
            uint w = (uint)input.GetLength(0);
            uint h = (uint)input.GetLength(1);
            uint i, j;
            double max, magnitude, log;

            double[,] logs = new double[w, h];
            double[,] output = new double[w, h];

            max = 0.0;
            for (i = 0; i < w; i++)
            {
                for (j = 0; j < h; j++)
                {
                    magnitude = input[i, j].Magnitude();
                    log = Math.Log(1 + magnitude);
                    if (log > max)
                        max = log;
                    logs[i, j] = log;
                }
            }

            for (i = 0; i < w; i++)
            {
                for (j = 0; j < h; j++)
                {
                    output[i, j] = logs[i, j] / max;
                }   
            }

            return output;
        }

        /// <summary>
        /// Obtient le module au carré d'une transformée de fourier 
        /// </summary>
        /// <param name="input">Table de <see cref="Complex"/> représentant la transformée de Fourier</param>
        /// <returns>Tableau de <see cref="Complex"/> représentant le module au carré de la transformée</returns>
        public static double[,] GetMagnitudeSquaredPlot(Complex[,] input)
        {
            uint w = (uint)input.GetLength(0);
            uint h = (uint)input.GetLength(1);
            uint i, j;
            double max, magnitude, log;

            double[,] logs = new double[w, h];
            double[,] output = new double[w, h];

            max = 0.0;
            for (i = 0; i < w; i++)
            {
                for (j = 0; j < h; j++)
                {
                    magnitude = input[i, j].MagnitudeSquared();
                    log = Math.Log(1 + magnitude);
                    if (log > max)
                        max = log;
                    logs[i, j] = log;
                }
            }

            for (i = 0; i < w; i++)
            {
                for (j = 0; j < h; j++)
                {
                    output[i, j] = logs[i, j] / max;
                }
            }

            return output;
        }

        /// <summary>
        /// Obtient la phase d'une transformée de fourier 
        /// </summary>
        /// <param name="input">Table de <see cref="Complex"/> représentant la transformée de Fourier</param>
        /// <returns>Tableau de <see cref="Complex"/> représentant la phase de la transformée</returns>
        public static double[,] GetPhasePlot(Complex[,] input)
        {
            uint w = (uint)input.GetLength(0);
            uint h = (uint)input.GetLength(1);
            uint i, j;
            double max, phase, log;

            double[,] logs = new double[w, h];
            double[,] output = new double[w, h];

            max = 0.0;
            for (i = 0; i < w; i++)
            {
                for (j = 0; j < h; j++)
                {
                    phase = input[i, j].Phase();
                    log = Math.Log(1 + phase);
                    if (log > max)
                        max = log;
                    logs[i, j] = log;
                }
            }

            for (i = 0; i < w; i++)
            {
                for (j = 0; j < h; j++)
                {
                    output[i, j] = logs[i, j] / max;
                }
            }

            return output;
        }
        #endregion

        #region Private Methods
        private static Complex[,] Compute(Complex[,] input, bool reverse)
        {
            Complex[,] output = input;
            Compute2D(ref output, reverse);
            return output;
        }

        private static void Compute2D(ref Complex[,] input, bool reverse)
        {
            uint w = (uint)input.GetLength(0);
            uint h = (uint)input.GetLength(1);
            uint i, j;
            Complex[] complex;

            if (reverse)
                Shift(ref input);
            
            complex = new Complex[w];
            for (j = 0; j < h; j++)
            {
                for (i = 0; i < w; i++)
                    complex[i] = input[i, j];

                Compute1D(ref complex, reverse);

                for (i = 0; i < w; i++)
                    input[i, j] = complex[i];
            }

            complex = new Complex[h];
            for (i = 0; i < w; i++)
            {
                for (j = 0; j < h; j++)
                    complex[j] = input[i, j];

                Compute1D(ref complex, reverse);

                for (j = 0; j < h; j++)
                    input[i, j] = complex[j];
            }

            if (!reverse)
                Shift(ref input);
        }

        private static void Compute1D(ref Complex[] input, bool reverse)
        {
            double delta, sine, angle;
            uint jump, step, group, pair, match, n;
            Complex product, factor, multiplier;

            ReverseBits(ref input);
            
            jump = 1;
            n = (uint)input.Length;
            angle = reverse ? 1.0 : -1.0;
            angle *= Math.PI;
            multiplier = new Complex(-1.0, 0.0);

            for (step = 1; step < n; step <<= 1)
            {
                jump = step << 1;
                delta = angle / (double)step;
                sine = Math.Sin(0.5 * delta);
                multiplier = new Complex(-2.0 * sine * sine, Math.Sin(delta));
                factor = new Complex(1.0, 0.0);
                for (group = 0; group < step; group++)
                {
                    for (pair = group; pair < n; pair += jump)
                    {
                        match = pair + step;
                        product = factor * input[match];
                        input[match] = input[pair] - product;
                        input[pair] += product;
                    }
                    factor = multiplier * factor + factor;
                }
            }

            if (!reverse)
                Scale(ref input);
        }

        private static void ReverseBits(ref Complex[] input)
        {
            uint i, j, k, i2;
            uint n = (uint)input.Length;
            Complex temp;

            i2 = n >> 1;
            j = 0;
            for (i = 0; i < n - 1; i++)
            {
                if (i < j)
                {
                    temp = input[i];
                    input[i] = input[j];
                    input[j] = temp;
                }
                k = i2;
                while (k <= j)
                {
                    j -= k;
                    k >>= 1;
                }
                j += k;
            }
        }

        private static void Scale(ref Complex[] input)
        {
            int n = input.Length;
            double factor = 1.0 / (double)n;
            for (int i = 0; i < n; i++)
            {
                input[i] *= factor;
            }
        }

        private static void Shift(ref Complex[,] input)
        {
            int i, j;
            int w = input.GetLength(0) / 2;
            int h = input.GetLength(1) / 2;
            Complex temp;

            for (i = 0; i < w; i++)
            {
                for (j = 0; j < h; j++)
                {
                    temp = input[i + w, j + h];
                    input[i + w, j + h] = input[i, j];
                    input[i, j] = temp;

                    temp = input[i + w, j];
                    input[i + w, j] = input[i, j + h];
                    input[i, j + h] = temp;
                }
            }
        }
        #endregion
    }
}
